% 猪肝样本的物理属性
rho = 1050; % 猪肝密度 (kg/m^3)
c = 3500; % 比热容 (J/kg·K)
kappa = 0.5; % 热导率 (W/m·K)

% 电极的加热区域体积 (椭球体积)
a = 0.01; % 长轴 (m)
b = 0.008; % 短轴 (m)
c_depth = 0.03; % 深度 (m)
V = (4/3) * pi * a * b * c_depth; % 椭球体积 (m^3)
P = 30; % 功率 (W)

% 离散化网格
r_max = 0.2; % 最大r
z_max = 0.2; % 最大z
dr = 0.001; % r的步长
dz = 0.001; % z的步长
dt = 0.01; % 时间步长
t_max = 8; % 最大时间

% 网格尺寸
r = 0:dr:r_max;
z = 0:dz:z_max;
t = 0:dt:t_max;

% 初始化温度场
T = zeros(length(r), length(z), length(t));

% 初始条件
T(:,:,1) = 25; % 初始温度 25℃

% 时间迭代计算
for n = 1:length(t)-1
    for i = 2:length(r)-1
        for j = 2:length(z)-1
            % 离散化热传导方程
            T(i,j,n+1) = T(i,j,n) + (kappa*dt/(rho*c)) * ((T(i+1,j,n) - 2*T(i,j,n) + T(i-1,j,n))/dr^2 + (T(i,j+1,n) - 2*T(i,j,n) + T(i,j-1,n))/dz^2) + (P/V) * dt;
        end
    end
end

time_step = 700;  % 选择展示的时间步
% 将不同的z位置上的温度场叠加
c = 3.587;
z_values = linspace(-c, c, 60);  % z的位置从-c到c，分成60个点
T_total = zeros(length(r), length(z));  % 初始化总温度场

% 叠加每个z位置的温度
for k = 1:length(z_values)
    z_pos = z_values(k);  % 当前z位置
    % 假设有一个函数来计算每个z位置的温度，下面是一个示例
    % 这里的calculate_temperature_at_z函数需要你根据实际情况进行编写
    T_at_z = calculate_temperature_at_z(r, z, z_pos);  % 计算当前位置的温度场
    T_total = T_total + T_at_z;  % 叠加温度场
end

% 绘制最终的温度场
figure;
surf(r, z, T_total);  % 绘制r-z平面上的温度分布
xlabel('r (m)');
ylabel('z (m)');
zlabel('Temperature (°C)');
title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
colorbar;  % 显示颜色条
axis equal;  % 保证坐标轴比例相同


% % 创建网格坐标
% % 创建网格坐标
% [R, Z] = meshgrid(r, z);
% 
% % 确保 T_at_time 是一个二维矩阵
% T_at_time = T(:,:,time_step);  % 提取指定时间点的温度数据
% 
% % 检查 T_at_time 的尺寸
% disp(size(T_at_time));  % 输出T_at_time的大小

% 使用 surf 绘制温度场
% figure;
% surf(R, Z, T_at_time);
% xlabel('r (m)');
% ylabel('z (m)');
% zlabel('Temperature (°C)');
% title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
% colorbar;  % 显示颜色条
% axis equal;  % 保证坐标轴比例相同
% figure;
% surf(R, Z, T_at_time);  % 绘制温度场
% shading interp;  % 平滑颜色过渡
% xlabel('r (m)');
% ylabel('z (m)');
% title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
% colorbar;  % 显示颜色条
% 
% % 移除 z 轴，隐藏 z 轴标签
% zlim([min(T_at_time(:)), max(T_at_time(:))]);  % 设置 z 轴的范围
% set(gca, 'ztick', []);  % 隐藏 z 轴刻度
% set(gca, 'zticklabel', []);  % 隐藏 z 轴刻度标签
% 
% axis equal;  % 保证坐标轴比例相同
% 定义r和z的区间和分辨率
% r_min = -0.2;
% r_max = 0.2;
% z_min = -0.2;
% z_max = 0.2;
% num_points = 100;  % 分成100格
% 
% % 创建r和z的网格
% r = linspace(r_min, r_max, num_points);  % 创建r的线性区间
% z = linspace(z_min, z_max, num_points);  % 创建z的线性区间
% [R, Z] = meshgrid(r, z);  % 生成r和z的网格

% 假设T_at_time是已知的温度数据矩阵，必须确保它与R和Z的大小匹配
% 这里假设T_at_time是大小为[100, 100]的二维矩阵
% T_at_time = ...  % 温度数据，你需要将它加载到这里

% 绘制温度分布的颜色映射
% 使用 imagesc 绘制颜色映射图
% figure;
% imagesc(r, z, T_at_time);  % 绘制二维温度场
% axis xy;  % 设置坐标轴方向正确
% xlabel('r (m)');
% ylabel('z (m)');
% title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
% colorbar;  % 显示颜色条
% 
% % 移除 z 轴，隐藏 z 轴标签
% zlim([min(T_at_time(:)), max(T_at_time(:))]);  % 设置 z 轴的范围
% set(gca, 'ztick', []);  % 隐藏 z 轴刻度
% set(gca, 'zticklabel', []);  % 隐藏 z 轴刻度标签
% 
% axis equal;  % 保证坐标轴比例相同
% 参数设置
% c = 3.587;         % 热源的半径（单位：m）
% num_points = 100;  % 每个方向的网格点数
% r_min = -0.2;      % r轴的最小值
% r_max = 0.2;       % r轴的最大值
% z_min = -c;        % z轴的最小值
% z_max = c;         % z轴的最大值
% 
% % 划分网格
% r = linspace(r_min, r_max, num_points);   % 在r方向上生成网格
% z = linspace(z_min, z_max, num_points);   % 在z方向上生成网格
% [R, Z] = meshgrid(r, z);  % 创建网格坐标
% 
% % 初始化最终的温度场
% T_final = zeros(num_points, num_points);
% 
% % 假设T_at_time是根据偏微分方程计算得到的每个点的温度分布
% % 可以根据你的热传导方程计算每个点的温度。以下是一个简化的假设。
% % 这里以T_at_time作为每个z位置下的温度分布矩阵，按位置叠加。
% 
% for k = 1:num_points
%     % 计算每个z点的温度场（这里只是一个简化的假设，具体根据偏微分方程计算）
%     % 你可以根据你的热传导方程来计算每个z位置的温度T_at_z
%     T_at_z = calculate_temperature_at_z(R, Z, z(k));  % 这里调用一个函数计算z位置的温度
%     T_final = T_final + T_at_z;  % 将每个z位置的温度叠加
% end
% 
% % 绘制最终温度场
% figure;
% surf(R, Z, T_final);
% shading interp;  % 平滑颜色过渡
% xlabel('r (m)');
% ylabel('z (m)');
% title('Temperature Distribution');
% colorbar;  % 显示颜色条
% axis equal;  % 保证坐标轴比例相同
% 
% 
